---
title: "Natural experiment"
author: "Kota SUECHIKA"
date: "`r format(Sys.time(), '%Y-%m-%d')`"
output: html_document
---

```{r setup, include=FALSE}
rm(list=ls())
knitr::opts_chunk$set(echo = TRUE, fig.width = 10, fig.height = 5)
require(memisc)
require(stringi)
library(stargazer)
library(coefplot)
library(interplot)
library(sjPlot)
library(sjmisc)
library(sjlabelled)
library(stargazer)
library(stringr)
library(tidyverse)
library(ggplot2)
library(psych)
library(dplyr)
library(modelsummary)
library(makedummies)
library(ggeffects)
library(estimatr)
library(ggpubr)
library(BalanceR)
library(scales)

theme_set(theme_sjplot())

dat_lebanon <- read_csv("data/Lebanon.csv")
dat_lebanon <- dat_lebanon[-1,]

dat_lebanon$RecordedDate <- as.Date(stri_match_first_regex(dat_lebanon$RecordedDate, "\\d{4}-\\d{2}-\\d{2}"))

dat_lebanon$ceasefire <- NA
dat_lebanon$ceasefire[dat_lebanon$RecordedDate <= as.Date("2025-01-16")] <- "Pre"
dat_lebanon$ceasefire[as.Date("2025-01-17") <= dat_lebanon$RecordedDate] <- "Post"
dat_lebanon$ceasefire <- factor(dat_lebanon$ceasefire)

dat_lebanon <- dat_lebanon |>
  mutate(Gender = D1,
         Age = D2,
         Education = D3,
         Income = D5,
         Hezbollah_support = Q4_3)

line_date <- as.Date("2025-01-16")
```

# timeseries sample distribution
```{r}
Sys.setlocale("LC_TIME", "en_US.UTF-8")

dailycount_lb <- dat_lebanon |>
  group_by(RecordedDate) |>
  summarise(Count = n()) |>
  arrange(RecordedDate)
dailycount_lb$RecordedDate <- as.Date(dailycount_lb$RecordedDate)

p1 <- ggplot(data = dailycount_lb, aes(x = RecordedDate, y = Count)) +
  geom_bar(stat = "identity", fill = "darkgrey") +
  geom_vline(xintercept = as.numeric(line_date), color = "black", linetype = "dashed", size = 1) +
  annotate("text", x = line_date, y = max(dailycount_lb$Count) * 1.05, 
           label = "Gaza Ceasefire", color = "black", angle = 90, vjust = -0.5, hjust = 1.2, size = 5) +
  scale_x_date(date_breaks = "5 days", date_labels = "%b %d") +
  labs(title = "Samples per day",
       x = "Date",
       y = "Number of samples") +
  theme(
    plot.title = element_text(size = 20, face = "bold"),  
    axis.title.x = element_text(size = 16),                
    axis.title.y = element_text(size = 16),                
    axis.text.x  = element_text(size = 14, angle = 45, hjust = 1), 
    axis.text.y  = element_text(size = 14)                  
  )

p1

ggsave(filename = "result/sample_per_day.png",
       plot = p1,
       width = 15,
       height = 7,
       units = "in", 
       dpi = 600,
       bg = "white")
```

# time series sample distribution refined
```{r}

p1_clean <- ggplot(data = dailycount_lb, aes(x = RecordedDate, y = Count)) +
  geom_bar(stat = "identity", fill = "gray40", width = 0.8) +
  geom_vline(xintercept = as.numeric(line_date), color = "black", linetype = "dashed", size = 0.8) +
  annotate("label", x = line_date, y = max(dailycount_lb$Count), 
           label = "Gaza Ceasefire Annoucement", fill = "white", color = "black", 
           vjust = 1, hjust = -0.1, size = 5) +
  scale_x_date(date_breaks = "1 week", date_labels = "%b %d") +
  
  labs(title = "",
       subtitle = "Number of respondents per day",
       x = "",
       y = "Count") +
  

  theme_classic(base_size = 16) +
  theme(
    plot.title = element_text(face = "bold", size = 20),
    axis.text.x = element_text(angle = 45, hjust = 1, color = "black"),
    axis.text.y = element_text(color = "black"),
    panel.grid.major.x = element_blank(),
    panel.grid.major.y = element_line(color = "gray90") 
  )

print(p1_clean)

ggsave("result/sample_per_day_clean.png", p1_clean, width = 12, height = 6, dpi = 600)
```

# balance check
```{r}
balance_lb <- BalanceR(data = dat_lebanon, 
                    group = ceasefire,
                    cov = c(Gender, Age, Education, Income))
print(balance_lb, digits = 3)
summary(balance_lb, digits = 5)
plot(balance_lb, point.size = 5, text.size = 18)

balance_LB <- plot(balance_lb, point.size = 5, text.size = 15) +
  labs(title = "Semi-natural experiment")
ggsave(filename = "result/balance_LB.png",
       plot = balance_LB,
       width = 7,
       height = 4,
       units = "in", 
       dpi = 600)
```

# regression
```{r}
reg10 <- lm(formula = Hezbollah_support ~ ceasefire, 
            data = dat_lebanon)
summary(reg10)

reg11 <- lm(formula = Hezbollah_support ~ ceasefire +
              Gender + Age + Education + Income, 
            data = dat_lebanon)
summary(reg11)

p_lebanon <- plot_model(reg10, type = "pred", terms = c("ceasefire"),
           axis.labels = "ceasefire", title = "Lebanon")
```

# plot
```{r}
ggsave(filename = "result/direct_question_result.pdf",
       plot = p_lebanon,
       width = 8,
       height = 5,
       units = "in", 
       dpi = 600)
```